Single-cell multiomics decodes regulatory programs for mouse secondary palate development

Perturbations in gene regulation during palatogenesis can lead to cleft palate, which is among the most common congenital birth defects. Here, we perform single-cell multiome sequencing and profile chromatin accessibility and gene expression simultaneously within the same cells (n = 36,154) isolated from mouse secondary palate across embryonic days (E) 12.5, E13.5, E14.0, and E14.5. We construct five trajectories representing continuous differentiation of cranial neural crest-derived multipotent cells into distinct lineages. By linking open chromatin signals to gene expression changes, we characterize the underlying lineage-determining transcription factors. In silico perturbation analysis identifies transcription factors SHOX2 and MEOX2 as important regulators of the development of the anterior and posterior palate, respectively. In conclusion, our study charts epigenetic and transcriptional dynamics in palatogenesis, serving as a valuable resource for further cleft palate research.

2. Page 4, line 29 and page 5, lines 1-2 refer to Figure S1C top and bottom, when the panels are presented as left and right.
3. The label "Glial" is missing from the y-axis in Figure 1E.4. The RNA and motif p-values for Twist2 on page 7 do not match the values in Table 1.

5.
It is unclear what the top and bottom panels represent at each location in Figure 3D.On page 8, lines 15-16 the authors state that "Inhba was restricted to beneath the epithelial layer in the anterior region", when a similar expression pattern is also noted in the middle region.Similarly, the authors state that "Sim2 and Trps1 were expressed only in the anterior half of the posterior region", though they appear to be expressed medially.
6. On page 8, the authors should define how the "anterior" and "posterior" regions of the secondary palate were delineated for bulk RNA-seq.
7. The authors should comment on why the differentially expressed genes with the largest fold changes in bulk RNA-seq in Figure 3F were often not detected in the scRNA-seq.Relatedly, it is unclear to this reviewer what is being represented in Figure 3H.Does this indicate that only ~6 transcripts commonly mapped to anterior and posterior locations between datasets?8.It is unclear why none of the labeled driver genes for posterior trajectory in Figure 5C appear in Figure 4L, especially Meox2.The colors in the legend of Figure 5C do not appear to match the colors in the graph.
Reviewer #3 (Remarks to the Author): The manuscript provides a comprehensive analysis of the developmental programs of the mouse secondary palate, utilizing single-cell multiome data to simultaneously profile RNA-seq and ATAC-seq in the same cells.The study identified distinct cell types and developmental trajectories, as well as key regulators and signaling pathways involved in the development of the mouse secondary palate.Overall, I find the paper to be a useful reference and data resource for studying developmental processes at single-cell resolution.
-The cell clustering and annotation appear to be driven solely by signals in RNA data, while recent computational methods (e.g.Seurat v4) can handle paired RNA and ATAC data to create joint embeddings of the two modalities for clustering.Therefore, using information from both modalities for clustering would be more reasonable.
-It would be helpful to understand how the authors selected the 3000 highly variable genes and whether the annotation results are stable to the chosen genes.
-Checking motif enrichment of ATAC data would be beneficial in addition to checking the enrichment of gene activity scores for differentially expressed genes in the annotated cell types.
-How consistent is the diffusion pseudotime against real time (days) in the trajectory analysis?-The driver genes were inferred by checking marginal correlation with the fate probability, but it should be acknowledged that not all marginally correlated genes are necessarily drivers.There could be many "passenger" genes correlated with the real driver genes, but counted as driver genes using this method.This limitation should be discussed.Response: We thank the reviewer for this valuable suggestion.We would like to point out that in our precious version, the CellRank algorithm used in Supplementary Figure 8 detects the initial and terminal cell states of the system and computes a global map of fate potentials, assigning each cell, including the progenitor populations, the probability of reaching each terminal state.
In addition, as suggested by the reviewer, we compared cell fate decisions within the progenitor populations directly.We added the description of these results to the main text and summarized in a new Supplementary Figure 9.In the revised manuscript, on page 10, we added: "To more granularly resolve these velocity predictions, we applied CellRank, which infers initial and terminal state probabilities for each cell based on RNA velocity.Consistent with WOT-derived trajectories, CellRank found high initial state probabilities in CNC-derived multipotent cells and high terminal state probabilities in the late-stage subpopulations (Supplementary Fig. 8).To analyze cell fate decisions within the progenitor populations, we isolated those with high transition probabilities toward the anterior and posterior trajectories, respectively.Differential gene expression analysis unveiled distinct expression profiles in these two subpopulations.In the progenitor anterior subpopulation, genes such as Shox2, Satb2, Inhba, Cyp26b1, and Nrp1 demonstrated significantly higher expression levels.Conversely, in the progenitor posterior subpopulation, genes like Meox2, Prickle1, Sim2, Efnb2, and Trps1 exhibited elevated expression (Supplementary Fig. 9).The results align with the expression profile observed in terminally differentiated anterior and posterior populations, as illustrated in Figure 3".Supplementary Fig. 9. Progenitor cells with high transition probabilities displayed higher anterior-specific genes expressions.a Scatter plot displays transition probabilities to anterior (xaxis) and posterior (y-axis) states for progenitor cells.b Volcano plot shows average log2 fold change (x-axis) and -log10 adjusted p-value (y-axis) of differentially expressed genes between progenitor cells biased towards the anterior and posterior cell fates.
3. Most importantly, the authors must use some approach to validate the most key findings in a functional experiment.For this, the authors shall use lentiviral microinjections, CAS9/CRISPR knockouts and analysis in F0 or F1 and so on, and also the regulatory regions reporter assays.Without extensively addressed point 3, this paper will likely fail to be published in Nature Communications.

Response:
We thank the reviewer for this critical point.We agree that functional assays are important.However, we were not able to perform CRISPR/Cas9 knockouts because generation of mouse strains by CRISPR/Cas 9 for the candidate genes/mutations will take many months, but we were given only three months.In addition, this is somewhat out scope of our current study, which has already been extensive.Please kindly note that, although single-cell multiome technologies are Supplementary Fig. 9. Progenitor cells with high transitio n probabilities displayed higher anterior-specific genes expressions.a Scatter plot of transition probability to anterior (x-axis) and posterior (y-axis) of progenitor cells.b Volcano plot shows average log2 fold change (x-axis) and -log10 adjusted p-value (y-axis) of dif ferentially expressed genes between progentor anterior cells and progenitor poster ior cells.Progenitor anterior vs Progentor posterior not very new, there has been no such a study in the field of craniofacial development, a critical development stage important for not only developmental biology but also human diseases.We sincerely hope that you will agree our explanation.Instead, to experimentally validate our key findings, we generated ChIP-seq experiments and included these in our revision.
On pages 11 and 12, we wrote "We next investigated underlying transcriptional regulators by motif enrichment analysis.We characterized a list of lineage-determining TFs that are potentially candidates to control each trajectory by binding to regulatory elements to regulate the expression of the above-mentioned driver genes.Shox2 was identified as an important regulator at the start of the anterior trajectory (motif adjusted p-value = 6.09x10-3, motif fold change = 4.21, gene adjusted pvalue < 2.2x10-16, gene fate correlation = 0.47) (Figure 5a).MEOX2 showed potential regulatory roles in the middle (motif adjusted p-value = 0.024, motif fold change = 2.87) and end of the posterior trajectory (motif adjusted p-value = 1.85x10-3, motif fold change = 2.06, Figure 5b).To experimentally validate these predictions, we conducted chromatin immunoprecipitation followed by sequencing (ChIP-seq) experiments.SHOX2 and MEOX2 ChIP-seq data were generated from anterior and posterior palate tissue, respectively.We calculated the likelihood of observing a ChIPseq peak near the genes that were upregulated at the start of the anterior trajectory and the middle of the posterior trajectory.As computationally predicted, the likelihood of observing a SHOX2 peak was increased for genes up-regulated at the start of the anterior trajectory.Correspondingly, the likelihood of observing a MEOX2 peak was increased for genes up-regulated at the middle of the posterior trajectory (Figure 5c).For example, we observed a ChIP-seq binding peak for SHOX2 but not MEOX2 in the predicted SHOX2 target Nfia (Figure 5d).Simultaneously, a MEOX2 binding peak was observed upstream of the predicted MEOX2 target Has2 (Figure 5e)."In addition, we re-analyzed published RNA-seq data of Shox2 knockout to validate our in silico perturbation predictions in the revision on page 13: "To validate the predicted gene-regulatory dynamics, we re-analyzed published bulk RNA-seq data derived from the E14.5 anterior hard palatal tissues of Shox2 Cre/-mice, wherein the Shox2 gene had been knocked out.Among the 11 predicted SHOX2 targets, 8 genes exhibited significantly altered expression following Shox2 knockout (adjusted p-value < 0.05) (Figure 6f).As expected, genes predicted to be positively regulated by Shox2 demonstrated decreased expression upon Shox2 knockout (Figure 6e, f).Correspondingly, genes predicted to be negatively regulated displayed increased expression upon Shox2 knockout (Figure 6e, f).Taken together these results suggest that SHOX2 and MEOX2 serve as crucial regulators driving the development of the anterior and posterior secondary palate, respectively."Together with several experiments included in the previous version (bulk RNA-seq, in situ hybridization, quantitative RT-PCR, etc.), we hope that the reviewer will agree that our main findings are validated and that this original work is a novel contribution in craniofacial development field.

Response to Reviewer #2
This manuscript details results obtained from combined single-cell RNA-seq and ATAC-seq of mouse secondary palatal shelves from E12.5-E14.5.The authors mostly focus on the cranial neural crest-derived mesenchyme, demonstrating differentiation into anterior palatal mesenchyme, posterior palatal mesenchyme, osteogenic, dental mesenchyme and perimysial lineages, consistent with previous findings.They further analyze the above data to predict transcription factors that determine individual lineages, such as Shox2 in the anterior palatal mesenchyme, again consistent with previous findings.The data is robust, clearly presented and will likely be a valuable resource for the field.However, it is unclear what new information is gained from this study, especially given the lack of validation studies and the fact that the majority of transcription factors discussed already have demonstrated in vivo roles in various secondary palate cell types.

Response:
We thank the reviewer #2 for carefully reading our manuscript and summarizing our work nicely.We understood some results of our work are consistent with previous findings, which demonstrated that our single-cell multiomics approach is reliable and informative.However, our study investigated craniofacial development at the single-cell multiomics level giving insight into the underlying mechanisms at unprecedented resolution.The extensive data and analysis will provide important knowledge and further understanding of mechanisms underlying craniofacial disease.
Due to the limited scope of the current work, we were not able to develop mouse strains or cell lines for additional experimental validation.Such experiments will take many months, but we were given only three months for revision.Please kindly note that, although single-cell multiome technologies are not very new, there has been no such a study in the field of craniofacial development, a critical development stage important for not only developmental biology but also human diseases.We sincerely hope that you will agree our explanation.Instead, we performed ChIP-seq experiments to validate our key findings (Shox2 and Meox2) and re-analyzed published RNA-seq of Shox2 knockout in the revision.Together with several experiments included in the previous version (bulk RNA-seq, in situ hybridization, quantitative RT-PCR, etc.), we hope that the reviewer will agree that our main findings are validated and that this original work is a novel contribution in craniofacial development field.
Major comments: 1. Page 4, line 20: Development of the mouse secondary palate occurs from E11.5-E16.5 (see Bush and Jiang, 2012, PMC3243091).Relatedly, the image of secondary palatal shelves fusing at E14.5 in Figure 1A is premature, as this representation more closely matches E15.0.The authors should clearly discuss why they chose the E12.5, E13.5, E14.0 and E14.5 timepoints and which steps of secondary palate development they were capturing at each timepoint.

Response:
We thank the reviewer for the valuable comment.We selected the specific timepoints based on the mouse anatomy described by the FaceBase project.Color-coded scanning electron microscopy images show the fusing of the palatal shelves at 14.5: Screenshot taken from FaceBase website (https://www.facebase.org/resources/mouse/mouseanatomy/)shows fusing of palatal shelves at E14.5.
Response: We apologize for not being clear.These markers were previously shown to be expressed in lateral or medial regions.In the manuscript by Lan et al., the authors stated "From E12.5 to E13.5, Osr2 mRNA is expressed abundantly throughout the palatal mesenchyme, with lateral regions expressing higher levels than the medial regions."Similarly, in the manuscript by Levi et al., the authors reported "Co-expression of Msx1 and Dlx5 is limited to the most lateral region of the maxillary mesenchyme."However, in the revised manuscript, we focused on anterior and posterior subpopulations, so the results related to medial and lateral positions were removed from the main text.4. The authors state on page 11 that "…we characterized a list of lineage-determining TFs that control the anterior trajectory, such as Shox2 at the early stage, Foxl2 at the middle stage, and Nr2e1 at the late stage of the trajectory", but do not provide any statistics on the latter two factors.Further, given the long list of transcription factors indicated at the right of Figure 4H, each with relatively low fold enrichment, it is hard to understand how the authors can make such a statement about Nr2e1, for example.To definitely back up such statements, the authors would need to perform ChIP-seq for transcription factors of interest and show that each binds at least a subset of the identified motifs, knock down individual transcription factor function in vivo and demonstrate that a specific trajectory is affected, and/or perform detailed lineage tracing experiments (for example, using novel Cre lines generated from this data).Without such data, which would significantly enhance the manuscript, the authors need to temper such statements as "…we identified a list of TFs that control each trajectory by binding to regulatory elements to regulate the expression of the above-mentioned driver genes" on page 14.

Response:
We thank the reviewer for this constructive suggestion.We agree that additional analysis was necessary to specify lineage determining transcription factors.In the revision, we added in silico perturbation analysis using CellOracle, a novel computational method to predict the impact of perturbation of transcription factors on cellular development.This analysis is included in the revised manuscript and the results are summarized in the new Figure 6 and Supplementary Figure 12 on pages 12 and 13: "In silico perturbation analysis analysis finds driver genes Next, we applied CellOracle to assess the impact of perturbing specific regulators on the development of the secondary palate.This algorithm leverages single-cell multi-omics data to deduce gene-regulatory networks.It then conducts in silico perturbations to simulate how these changes affect cellular development, relying solely on unperturbed data.
Our analysis focused on the cells within the anterior and posterior trajectories.Independent data analysis was conducted, including normalization, clustering, and dimension reduction using PAGA and force-directed graphs, followed by diffusion pseudotime calculation.The manifold revealed two distinct trajectories originating from the multipotent cells towards the anterior and posterior cells (Figure 6a).We then calculated perturbation scores for all detected TFs. High perturbation scores indicate that in silico knockout of the TF significantly decreased the development of the trajectory, suggesting that the TF is an essential driver of the trajectory.Interestingly, while the CellOracle perturbation scores were correlated for many TFs, SHOX2 and MEOX2 showed relatively high specificity for the anterior and posterior trajectories, respectively (Figure 6b).Indeed, in silico perturbation of Shox2 and Meox2 reversed the developmental velocities for the anterior and posterior trajectories, respectively (Figure 6c, d, Supplementary Fig. 13).
The regulatory networks predicted by CellOracle for all transcription factors can be found in Supplementary Table 3.For SHOX2 and MEOX2, CellOracle predicted 11 and 4 target genes, respectively.It is noteworthy that among these predicted targets, Satb2, Prrx1, and Prickle1 have previously been linked to cleft palate (Figure 6e)."In addition, we re-analyzed published RNA-seq data of Shox2 knockout to validate our in silico perturbation predictions in the revision on page 13: "To validate the predicted gene-regulatory dynamics, we re-analyzed published bulk RNA-seq data derived from the E14.5 anterior hard palatal tissues of Shox2 Cre/-mice, wherein the Shox2 gene had been knocked out.Among the 11 predicted SHOX2 targets, 8 genes exhibited significantly altered expression following Shox2 knockout (adjusted p-value < 0.05) (Figure 6f).As expected, genes predicted to be positively regulated by Shox2 demonstrated decreased expression upon Shox2 knockout (Figure 6e, f).Correspondingly, genes predicted to be negatively regulated displayed increased expression upon Shox2 knockout (Figure 6e, f).Taken together these results suggest that SHOX2 and MEOX2 serve as crucial regulators driving the development of the anterior and posterior secondary palate, respectively."Finally, we also performed ChIP-seq experiments, as suggested by the reviewer, to validate our key findings (Shox2 and Meox2).On page 12, we wrote "To experimentally validate these predictions, we conducted chromatin immunoprecipitation followed by sequencing (ChIP-seq) experiments.SHOX2 and MEOX2 ChIP-seq data were generated from anterior and posterior palate tissue, respectively.We calculated the likelihood of observing a ChIP-seq peak near the genes that were upregulated at the start of the anterior trajectory and the middle of the posterior trajectory.As computationally predicted, the likelihood of observing a SHOX2 peak was increased for genes upregulated at the start of the anterior trajectory.Correspondingly, the likelihood of observing a MEOX2 peak was increased for genes up-regulated at the middle of the posterior trajectory (Figure 5c).For example, we observed a ChIP-seq binding peak for SHOX2 but not MEOX2 in the predicted SHOX2 target Nfia (Figure 5d).Simultaneously, a MEOX2 binding peak was observed upstream of the predicted Meox2 target Has2 (Figure 5e)."  anterior region and Inhba was restricted to beneath the epithelial layer in the anterior and middle region.In contrast, Meox2, Prickle1, and Efnb2 were mainly expressed in the posterior region of the palate.Interestingly, Sim2 and Trps1 were expressed medially in the anterior half of the posterior region of the developing secondary palate." 6. On page 8, the authors should define how the "anterior" and "posterior" regions of the secondary palate were delineated for bulk RNA-seq.

Response:
We thank the reviewer for pointing this out.In the revised manuscript, we have addressed this by providing a clearer explanation of how we delineated the "anterior" and "posterior" regions of the secondary palate for bulk RNA-seq analysis.
In the revised manuscript, on page 8, we have updated the relevant sentence "To validate the gene signatures of anterior and posterior subpopulations, we performed bulk RNA-sequencing (RNA-seq) after isolating RNA from the microdissected anterior 1/3 (n=3) and posterior 1/3 regions (n=3) of the developing secondary palate at E14.0".
On page 18, we made the following revision: "The anterior (n=3) and posterior (n=3) 1/3 palatal shelves of E14.0 C57BL/6J mice were microdissected, isolated, and then subjected to bulk RNA sequencing." 7. The authors should comment on why the differentially expressed genes with the largest fold changes in bulk RNA-seq in Figure 3F were often not detected in the scRNA-seq.Relatedly, it is unclear to this reviewer what is being represented in Figure 3H.Does this indicate that only ~6 transcripts commonly mapped to anterior and posterior locations between datasets?
Response: We thank the reviewer for pointing this out.We examined the differentially expressed genes with large fold changes in bulk RNA-seq but were not detected in the scRNA-seq.We found the sequencing depth of these genes is significantly lower than the random gene set of the same length (p = 9.4e-13).In the revised manuscript, we added a new Supplementary Figure 6 and on page 8, we wrote "Although the highly significant correspondence between bulk and single-cell levels, several differentially expressed genes with large fold changes in bulk RNA-seq were not detected in the scRNA-seq data.These genes tended to be expressed at low levels in the bulk RNAseq, suggesting that the limited sensitivity of lowly expressed transcripts in the scRNA-seq assay may obscure the signals (Supplementary Fig. 6)".
Response: We thank the reviewer for this important question.In the revised version, we added a new Supplementary figure panel depicting the relation between diffusion pseudotime and real time.
The increasing trend over time shows the consistency between these two variables.
In the revised manuscript, on page 10, we added "The cells with large diffusion pseudotime values tended to be derived from late time points (Supplementary Fig. 10)." Our analysis focused on the cells within the anterior and posterior trajectories.Independent data analysis was conducted, including normalization, clustering, and dimension reduction using PAGA and force-directed graphs, followed by diffusion pseudotime calculation.The manifold revealed two distinct trajectories originating from the multipotent cells towards the anterior and posterior cells (Figure 6a).We then calculated perturbation scores for all detected TFs. High perturbation scores indicate that in silico knockout of the TF significantly decreased the development of the trajectory, suggesting that the TF is an essential driver of the trajectory.Interestingly, while the CellOracle perturbation scores were correlated for many TFs, SHOX2 and MEOX2 showed relatively high specificity for the anterior and posterior trajectories, respectively (Figure 6b).Indeed, in silico perturbation of Shox2 and Meox2 reversed the developmental velocities for the anterior and posterior trajectories, respectively (Figure 6c, d, Supplementary Fig. 13).
The regulatory networks predicted by CellOracle for all transcription factors can be found in Supplementary Table 3.For SHOX2 and MEOX2, CellOracle predicted 11 and 4 target genes, respectively.It is noteworthy that among these predicted targets, Satb2, Prrx1, and Prickle1 have previously been linked to cleft palate (Figure 6e).
To validate the predicted gene-regulatory dynamics, we re-analyzed published bulk RNA-seq data derived from the E14.5 anterior hard palatal tissues of Shox2 Cre/-mice, wherein the Shox2 gene had been knocked out.Among the 11 predicted Shox2 targets, 8 genes exhibited significantly altered expression following Shox2 knockout (adjusted p-value < 0.05) (Figure 6f).As expected, genes predicted to be positively regulated by Shox2 demonstrated decreased expression upon Shox2 knockout (Figure 6e, f).Correspondingly, genes predicted to be negatively regulated displayed increased expression upon Shox2 knockout (Figure 6e, f).Taken together these results suggest that SHOX2 and MEOX2 serve as crucial regulators driving the development of the anterior and posterior secondary palate, respectively."Response: We thank the reviewer #2 for carefully evaluating our revision and the valuable advice.While the importance of Shox2 and Meox2 in secondary palate development has been demonstrated before, the precise lineage-specific regulatory networks at single-cell resolution have not yet been described before.In addition, while our manuscript focused on these two exemplary transcription factors, our data contains the predicted regulatory networks for a much larger collection of transcription factors representing a unique resource to the research community.We highlighted these two transcription factors to illustrate the useful information of our data.In the Results section, on page 16, we wrote: "While our analysis focused on these two important regulators of secondary palate development, our results provide additional information at genome-scale using single cell multiome approach.In addition to Shox2 and Meox2, we provide inferred regulatory networks for several additional transcription factors, which represents a valuable resource to the research community (Supplementary Table 3)." 1.The authors still have not discussed which steps of secondary palate development they were capturing at each timepoint.This will be crucial to make the findings clear to a broad audience interested in tissue morphogenesis.To be clear, the authors should be discussing milestones such as downgrowth, elevation, fusion, etc.

Response:
We thank the reviewer for this critical point and constructive suggestion.We apologize for any lack of clarity in the previous revision.To clarify which steps of secondary palate development were captured at each timepoint, we added the following text to the second revised manuscript.
In the Methods section, on page 17, we wrote: "Palatal shelves were isolated from time-mated C57BL/6J mice (000664, Jackson Laboratory) at four distinct developmental timepoints, including E12.5, E13.5, E14.0, and E14.5.Specifically, at E12.5, we primarily captured the early stages of palatal shelves, during which they emerge as outgrowths from the maxillary processes.E13.5 corresponds to a phase characterized by palatal shelf downgrowth towards the tongue.At E14.0, our samples included palatal shelves in the process of Supplemental Fig. 15.SCENIC+ identified regulons for anterior and posterior subpopulations, respectively.Scatter plot shows regulon specificity score (RSS, y-axis) of ranked regulons (x-axis) in anterior (left) and posterior subpopulations (right).
Regarding the Meox2 regulatory network, in the previous version, we employed a stringent p-value threshold (p-value < 0.001) to filter the TF-target links, resulting in a relatively small predicted Meox2 regulatory network.In the revised manuscript, we included a larger regulatory network resulting from less stringent, but still significant, thresholding (p-value < 0.01).This expansion enabled the identification of additional Meox2 targets.Importantly, a significant fraction of these targets was previously associated with cleft palate (p-value = 1.01x10 -4 ).The following text was added to the Results section on page 14: "We expanded the MEOX2 network, by utilizing a less stringent p-value threshold (p-value < 0.01) and identified a total of 39 target genes.Next, we integrated this list of genes with the CleftGeneDB database 47 which collects curated genes with experimental evidence for relevance in cleft palate.Importantly, 6 of these 39 MEOX2 targets have previously been associated with cleft palate, which represents a significantly larger overlap than expected by chance (Fisher's Exact test, p-value = 0.0002, odds ratio = 9.2, Supplementary Fig. 16a).For instance, Cacna1d has been reported to be related to oral cleft phenotype in the GWASdb SNP-Phenotype Associations dataset 48 .Additionally, Foxp2 is linked to nonsyndromic cleft lip and/or palate through genome-wide linkage analysis 49 .Notably, our MEOX2 CHIP-seq data provided further evidence of MEOX2 binding to promoter regions of these genes (Supplementary Fig. 16b)."Following the reviewer's suggestion, we performed additional analysis to quantify the statistical significance of the enrichment between the scATAC peaks and the CHIP-seq data.On page 12, we wrote: "Despite limited statistical power given the use of duplicates, these results showed marginal significance (One-sided t-test, anterior start p-value = 0.055, posterior middle p-value = 0.09).For example, we observed a ChIP-seq binding peak for SHOX2 but not MEOX2 in the predicted SHOX2 target Nfia (p-value = 2.29x10 -4 , signal = 3.01, Figure 5d).Simultaneously, a MEOX2 binding peak was observed upstream of the predicted MEOX2 target Has2 (p-value = 1.53x10 -4 , signal = 2.82, Figure 5e)." The Figure 5c was also updated.In addition, we have adjusted the language in the manuscript to temper statements about the roles of Shox2 and Meox2, referring to them as important regulators rather than drivers.
The major innovation of our study is the derivation of high confidence gene regulatory networks important for the development of mouse secondary palate.We have validated our predicted TFtarget relationships through (1) newly generated ChIP-seq experiments of SHOX2 and MEOX2 binding as well as (2) published Shox2-knockout RNA-seq and published H3K27 acetylation experiments.While previous studies 1,2,3,4 have described the relevance of Shox2 and Meox2 in secondary palate biology, to the best of our knowledge our study is the first to derive highconfidence regulatory networks.While our manuscript primarily focuses on Shox2 and Meox2, our dataset includes results for several other genes, making it a valuable resource to the community.We understand the reviewer's point regarding the in vivo knockdown experiments, and we agree that such experiments are valuable.However, the scope of the current work limited our ability to conduct these experiments within the timeframe of this study.We will extend this work for future study which will have specific focus on the regulation mechanisms of these transcriptional factors in palate development.We hope that the additional validations, data integration, and analyses we have performed address the reviewer's concerns and demonstrate the significance of our findings in secondary palate development.
We sincerely appreciate the constructive feedback from Reviewer #2, which has greatly contributed to the refinement and clarity of our manuscript.We remain committed to providing robust and valuable insights into the regulatory programs governing secondary palate development.
In the Acknowledgement section, we have acknowledged the valuable suggestions from all three reviewers who helped improve the manuscript.

Figure 6 .
Figure 6.In silico perturbation analysis reveals SHOX2 and MEOX2 as drivers of the anterior and posterior trajectories, respectively.a Left: UMAP visualization of CNC-derived mesenchymal subpopulations.Right: Force directed graph visualizations of anterior and posterior subpopulations.b Scatter plot shows in silico TF perturbation scores for the anterior (x-axis) and posterior (y-axis) trajectories, respectively.c Left: CellOracle vector field graphic shows the developmental flow of anterior trajectory.Arrows start from multipotent cells and point towards the anterior subpopulation.Right: CellOracle vector field graphic shows simulated vector shift after in silico Shox2 knockout.The developmental flow towards anterior cells is reversed upon in silico knockout of Shox2.d Same visualization as in panel c, showing simulated Meox2 perturbation in the posterior trajectory.e Network graphs visualized using Cytoscape show predicted SHOX2 and MEOX2 regulatory networks.The shape of the arrow heads is based on the predicted direction of regulation.Target nodes are colored based on regulation in Shox2-knockout RNA-seq data.f Heatmap shows log2 fold change of RNA expression in three Shox2 knockout samples normalized to controls.The bulk RNAseq data used here was downloaded from Xu et al. 2019.

Figure 5 .
Figure 5. Key transcription factors driving cells towards anterior and posterior states are validated in ChIP-seq experiments.a, b Dot plots show enriched motifs (y-axis) at different stages of the (a) anterior trajectory (x-axis) and (b) posterior trajectories.Dots are scaled by motif enrichment ratio and colored by significance.cLeft, boxplot shows the likelihood of peaks mapping near genes upregulated at the start of the anterior trajectory (y-axis) for SHOX2 and MEOX2 binding (x-axis).Right, boxplot shows the likelihood of peaks mapping near genes upregulated at the middle of the posterior trajectory (y-axis) for SHOX2 and MEOX2 binding (x-axis).d Integrative Genomics Viewer (IGV) view depicts ChIP-seq binding peak for SHOX2 in the predicted SHOX2 target Nfia.e IGV view presents ChIP-seq binding peak for MEOX2 upstream of the predicted MEOX2 target Has2.

Figure 6 .
Figure 6.In silico perturbation analysis reveals SHOX2 and MEOX2 as drivers of the anterior and posterior trajectories, respectively.…f Heatmap shows log2 fold change of RNA expression for predicted targets in three Shox2 knockout samples normalized to controls.The bulk RNA-seq data used here was downloaded from Xu et al. 2019.

Figure 1 .
Figure 1.Single-cell multiome assays dissect transcriptome and epigenome changes of the developing mouse secondary palate.…(d, e) Dot plot illustrates marker gene expression (x-axis) (d) and gene activity (e) (x-axis) across cell types (y-axis).Dot size is proportional to the percent of expressed cells.Colors indicate low (purple) to high (yellow) expression.

Figure 6 .
Figure 6.In silico perturbation analysis reveals SHOX2 and MEOX2 as drivers of the anterior and posterior trajectories, respectively.a Left: UMAP visualization of CNC-derived mesenchymal subpopulations.Right: Force directed graph visualizations of anterior and posterior subpopulations.b Scatter plot shows in silico TF perturbation scores for the anterior (x-axis) and posterior (y-axis) trajectories, respectively.c Left: CellOracle vector field graphic shows the developmental flow of anterior trajectory.Arrows start from multipotent cells and point towards the anterior subpopulation.Right: CellOracle vector field graphic shows simulated vector shift after in silico Shox2 knockout.The developmental flow towards anterior cells is reversed upon in silico knockout of Shox2.d Same visualization as in panel c, showing simulated Meox2 perturbation in the posterior trajectory.e Network graphs visualized using Cytoscape show predicted SHOX2 and MEOX2 regulatory networks.The shape of the arrow heads is based on the predicted direction of regulation.Target nodes are colored based on regulation in Shox2-knockout RNA-seq data.f Heatmap shows log2 fold change of RNA expression in three Shox2 knockout samples normalized to controls.The bulk RNAseq data used here was downloaded from Xu et al. 2019.

Figure 6 .
Figure 6.In silico perturbation analysis reveals SHOX2 and MEOX2 as drivers of the anterior and posterior trajectories, respectively.…f Heatmap shows log2 fold change of RNA expression for predicted targets in three Shox2 knockout samples normalized to controls.The bulk RNA-seq data used here was downloaded from Xu et al. 2019.

Figure 5 .
Figure 5. Key transcription factors driving cells towards anterior and posterior states are validated in ChIP-seq experiments.a, b Dot plots show enriched motifs (y-axis) at different stages of the (a) anterior trajectory (x-axis) and (b) posterior trajectories.Dots are scaled by motif enrichment ratio and colored by significance.cLeft, boxplot shows the likelihood of peaks mapping near genes upregulated at the start of the anterior trajectory (y-axis) for SHOX2 and MEOX2 binding (x-axis).Right, boxplot shows the likelihood of peaks mapping near genes upregulated at the middle of the posterior trajectory (y-axis) for SHOX2 and MEOX2 binding (x-axis).d Integrative Genomics Viewer (IGV) view depicts ChIP-seq binding peak for SHOX2 in the predicted SHOX2 target Nfia.e IGV view presents ChIP-seq binding peak for MEOX2 upstream of the predicted MEOX2 target Has2.
11.The H3K27 acetylation profiles validate the inferred anterior and posterior trajectories.Bar plot shows overlap (y-axis) between the scATAC anterior and posterior peaks and the corresponding anterior and posterior acetylation tracks (x-axis).

Figure 5 .
Figure 5. Key transcription factors driving cells towards anterior and posterior states are validated in ChIP-seq experiments.…c Left, boxplot shows the likelihood of peaks mapping near genes upregulated at the start of the anterior trajectory (y-axis) for MEOX2 and SHOX2 binding (xaxis).Right, boxplot shows the likelihood of peaks mapping near genes upregulated at the middle of the posterior trajectory (y-axis) for MEOX2 and SHOX2 binding (x-axis).